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Abstract. Polynomial approximations of computationally intensive models are central to uncertainty quan- 
tification. This paper describes an adaptive method for non-intrusive pseudospectral approximation, based on 
Smolyak's algorithm with generalized sparse grids. We rigorously develop and generalize the non-adaptive method 
proposed in [6] , and compare it to a common alternative approach for using sparse grids to construct polynomial 
approximations, direct quadrature. Analysis of direct quadrature shows that C(l) errors are an intrinsic property 
of some configurations of the method, as a consequence of internal aliasing. We provide precise conditions, based 
on the chosen polynomial basis and quadrature rules, under which this aliasing error occurs. We then estab- 
lish theoretical results on the accuracy of Smolyak pseudospectral approximation, and show that the Smolyak 
approximation avoids internal aliasing and makes far more effective use of sparse function evaluations. These 
results are applicable to broad choices of quadrature rule and generalized sparse grids. Exploiting this flexibility, 
we introduce a greedy heuristic for adaptive refinement of the pseudospectral approximation. We numerically 
demonstrate convergence of the algorithm on the Genz test functions, and illustrate the accuracy and efficiency 
of the adaptive approach on a realistic chemical kinetics problem. 
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1. Introduction. A central issue in the field of uncertainty quantification is the response 
of a model to random inputs. When model evaluations are computationally intensive, techniques 
for approximating the model response in an efficient manner are essential. Approximations may 
be used to evaluate moments or the probability distribution of a model's outputs, or to evaluate 
sensitivities of model outputs with respect to the inputs [20, 41, 33]. Approximations may also 
be viewed as surrogate models to be used in optimization [21] or inference [22], replacing the full 
model entirely. 

Often one is faced with black box models that can only be evaluated at designated input 
points. We will focus on constructing multivariate polynomial approximations of the input- 
output relationship generated by such a model; these approximations offer fast convergence for 
smooth functions and are widely used. One common strategy for constructing a polynomial ap- 
proximation is interpolation, where interpolants are conveniently represented in Lagrange form 
[1, 42]. Another strategy is projection, particularly orthogonal projection with respect to some 
inner product. The results of such a projection are conveniently represented with the correspond- 
ing family of orthogonal polynomials [4, 20, 43]. When the inner product is chosen according to 
the input probability measure, this construction is known as the (finite dimensional) polynomial 
chaos expansion (PCE) [14, 32, 8]. Interpolation and projection are closely linked, particularly 
when projection is computed via discrete model evaluations. Moreover, one can always realize 
a change of basis [9] for the polynomial resulting from either operation. Here we will favor or- 
thogonal polynomial representations, as they are easy to manipulate and their coefficients have 
a useful interpretation in probabilistic settings. 

This paper discusses adaptive Smolyak pseudospectral approximation, an accurate and com- 
putationally efficient approach to constructing multivariate polynomial chaos expansions. Pseu- 
dospectral methods allow the construction of polynomial approximations from point evaluations 
of a function [4, 3]. We combine these methods with Smolyak's algorithm, a general strategy 
for sparse approximation of linear operators on tensor product spaces, which saves computa- 
tional effort by weakening the assumed coupling between the input dimensions. Gerstner and 
Griebel [13] developed an adaptive variant of Smolyak's algorithm for numerical integration and 
illustrated the effectiveness of on-the-fly heuristic adaptation. Here we extend their approach 
to the pseudospectral approximation of functions. Adaptivity is expected to yield substantial 
efficiency gains in high dimensions — particularly for functions with anisotropic dependence on 
input parameters and functions whose inputs might not be strongly coupled at high order. 
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Previous attempts to extend pseudospectral methods to multivariate polynomial approxima- 
tion with sparse model evaluations employed ad hoc approaches that are not always accurate. 
A common procedure has been to use sparse quadrature, or even dimension-adaptive sparse 
quadrature, to evaluate polynomial coefficients directly [39, 20]. This leads to at least two dif- 
ficulties. First, the truncation of the polynomial expansion must be specified independently of 
the quadrature grid, yet it is unclear how to do this, particularly for anisotropic and generalized 
sparse grids. Second, unless one uses excessively high-order quadrature, significant aliasing errors 
may result. Constantine et al. provided the first clear demonstration of these aliasing errors, and 
proposed a Smolyak algorithm that does not share them [6]. 

Our treatment places this solution in the broader context of Smolyak constructions, and 
explains the origin of the observed aliasing errors for general (e.g., anisotropic) choices of sparse 
grid and quadrature rule. The main contributions of this paper are twofold. The first is theo- 
retical: we provide an error analysis that shows why the conventional approach (termed 'direct 
quadrature') goes wrong, and we develop polynomial exactness results for Smolyak pseudospec- 
tral approximation on generalized sparse grids. We do so without relying on properties of any 
particular quadrature rule, beyond univariate polynomial exactness. These results then provide a 
rigorous foundation for adaptivity: we introduce and demonstrate a fully adaptive algorithm for 
Smolyak pseudospectral approximation, which uses a single tolerance parameter to drive iterative 
refinement of both the polynomial approximation space and the collection of model evaluation 
points. 

The remainder of this introduction provides a brief motivating example. Section 2 defines 
the essential building-block problems of integration and polynomial approximation. Section 3 
develops computable one-dimensional and tensorized approximations for these problems, i.e., 
quadrature and pseudospectral approximation. Section 4 describes general Smolyak algorithms 
and their properties, then applies these results to quadrature and pseudospectral approxima- 
tion, yielding our principal theorems about exactness of the Smolyak approximation. Section 5 
compares the Smolyak approach to conventional direct quadrature. We show that, in almost all 
cases, direct quadrature is not an appropriate method for constructing polynomial expansions 
and should be superseded by Smolyak pseudospectral methods. In Section 6 we present an adap- 
tive extension to the Smolyak approach. Sections 7 and 8 contain observations on the choice of 
quadrature rule and connections with polynomial interpolation. As the adaptive methods are 
heuristic. Section 9 provides numerical demonstrations of their usefulness. 

1.1. Overview. Before a more formal exposition, we give a brief overview of the Smolyak 
approach and the results of our analysis. This section is intentionally loose with notation, pre- 
ferring to focus on intuition. 

A polynomial chaos expansion is a series expansion: 



in some basis of orthonormal polynomials ^'^ [38, 43]. If the function / is square integrable, the 
series converges to / in mean square. Practical uses of polynomial chaos select some truncation 
of the series. Taking the inner product of both sides with a basis function yields an expression 
for the coefficients 



where Sij is the Kronecker delta. The pseudospectral approach approximates this expression by 
substituting a quadrature operator Q for the inner product [4] : 




(1.1) 




(1.2) 




(1.3) 
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This approach incurs aHasing errors because, while the basis functions are analytically orthogonal, 
some pairs of basis functions will not appear orthogonal when the inner product is approximated 
via quadrature. Specifically, if for some Q(^'i^'j) ^ Sij, we say that ^'i and ^'j are not 
numerically orthogonal. These quadrature inaccuracies cause aliasing, because one coefficient 
can "bleed" onto another. 

Example 1.1. Assume that f{x) = ^2{x) + 10^^^'2o(x), so that /2 = 1 and /20 = 10^^. 
For the sake of this example, define three quadrature inner products as: 

Q(*2*2) = 1 (1.4) 

Q(«'20*20) = 1 (1.5) 
Q(*2*20) = 0.1 (1.6) 

Hence, the quadrature computes the mixed inner product incorrectly. If we compute the coeffi- 
cients with this quadrature rule, we find 

/2 = 1 * 1 + 0.01 * 10-3 = 1.0001 (1.7) 
/20 = 10-3 * 1 + 0.1 * 1 = .101 (1.8) 

Although /2 is essentially correct, /20 is not even the correct order of magnitude; hence this 
approximation is likely unacceptable. 

This example illustrates that it is crucial for pseudospectral methods to use quadrature that 
respects the orthogonality of the basis. In general, aliasing is an unavoidable consequence of 
the numerical approximation to the integral, so our goal is to limit its impact on the overall 
approximation. In one dimension, controlling the error is simple, but in higher dimensions the 
problem is more subtle. The nai've approach, direct quadrature, simply suggests using a sparse 
quadrature rule in (1.3). This approach can lose numerical orthogonality of the output basis, 
with disastrous results. This is precisely the source of error encountered by Constantine et al. 
[6]. 

The Smolyak algorithm resolves the potential for large errors by blending different full tensor 
approximations, which themselves are rigorously constructed to avoid certain aliasing errors; thus 
the Smolyak approximation inherits their favorable properties. The result is an algorithm that 
uses a sparse set of model evaluations, yet avoids the aliasing issues incurred by sparse direct 
quadrature schemes. 

2. One Dimensional and Tensor Problems. We begin by presenting problems of in- 
tegration and of approximation with polynomial chaos expansions. Exploiting the structure of 
tensor product problems is a core theme of this paper, so we will start in each case with a one- 
dimensional problem and then explain how to construct the corresponding multi-dimensional 
problem through tensorization. This section focuses on defining the problems and setting nota- 
tion, leaving numerical approximations to the next section. 

2.1. General setting. Consider a collection of one-dimensional linear operators where 
(i) indexes the operators used in different dimensions.^ We can extend a collection of these 
operators into higher dimensions by constructing the tensor product operator 

C'^'^'^ -.^ C^-^^ (g) ■ ■ ■ (g, C^^l (2.1) 

The one-dimensional operators need not be identical; the properties of the resulting tensor op- 
erator are constructed independently from each dimension. The bold parenthetical superscript 
refers to the tensor operator instead of the constituent one-dimensional operators. 



^By 'one-dimensional' we mean linear operators on functions of one variable, i.e., acts on functions of a;^*^. 
The tensor product operator in (2.1) then acts on functions of d variables, i.e., functions of x = {x^^\ . . . , x''*'). 
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2.2. One-dimensional integration. Let X^"^^ be an open or closed interval of the real line 
Then we define the weighted integral operator in one dimension as follows: 



dx 



where / : X(') 
function. 



is some real- valued function and w*^*-* : X^'^ 



(2.2) 

l^' is an integrable weight 



2.3. Multi-dimensional integration. One-dimensional integration extends easily to higher- 
dimensional integrals over product domains with separable weig ht functions. Let X^^\ X^"^) 
be a collection of intervals on the real line, as in the one dimensional case. Then let 



X := 



X(d) c 



(2.3) 



be the domain of integration, defined by the Cartesian product of one-dimensional intervals. 
Consider a real- valued function defined on this product space, / : X — >■ M. Then the multi- 
dimensional integral is given by the tensor product of one-dimensional integral operators: 

= 1(1) , 



.®IW(/) 



. w'^'^\x^'^^)f{x^^\. . . , x'^'^^) . . . dx^^^ (2.4) 



X 



where ?«(x) :— u)'-*-'(x'-'''). 



i=l 



2.4. One-dimensional polynomial chaos expansions. A polynomial chaos expansion 
approximates a function with a weighted sum of orthonormal polynomials [38, 43]. Let "HS^^ := 
{/ : ^ R} be a separable Hilbert space of square-integrable functions, with inner product 
for /, (? € H*^*^ defined as: 



(/,5> 



f{x)g{x)w'^'\x)dx^l'^'\fg). 



(2.5) 



Further, let the weight function be normalized such that (1, 1) = 1, so that w^'^'{x) may represent 
a probability density. There exist a set of polynomials {%Ij^^\x) : j e No} orthonormal with 
respect to this weight. Here No is the set of natural numbers including zero, and orthonormality 



requires {ip'j'\x),'ipk' (x)) — Sjk, yj,k e No. Without loss of generality, we restrict our attention 
to normalized polynomials. 

Let P„ be the space of univariate polynomials of degi 
an orthogonal projector onto this subspace. Then 

n 

-Pi'Hf) ■.= J2{fix),^f{x))i:f\x) 



j=0 



2=0 



^ pi*' be 



(2.6) 



with coefficients fj defined according to the middle term above. The orthogonal projection 
operator yields the optimal approximation of / in Pi*'', since any error is orthogonal to P^^ 
The polynomials are dense in H^'-* = (^X^^\w^^'>) , so the polynomial expansion of any f ^'H'^^^ 
converges in the sense as n — oo [4, 43]. If / S the coefficients must satisfy YlTLo ff < 

Projections with finite degree n omit terms of the infinite series, thus incurring truncation 
error. We can write this error as 



oo 



3=n+l 



OO 



j=n+l 



(2.7) 
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Hence, we may reduce the truncation error to any desired level by increasing n, removing terms 
from the sum in (2.7) [4, 16]. 

2.5. Tensor polynomial chaos expansions. Let 'H'-'^-' be the tensor product of the Hilbert 
spaces defined above: 

H^d) (2.8) 

'H^'^^ is a Hilbert space of real-valued functions on X, with inner product given by the integral 
form in (2.4), (/,<?) := l'^'^\fg). A basis {^'i : i G Nq} for this space may be formed by taking 
the tensor product of univariate orthonormal polynomials in each dimension, such that 

d 

vl/i(x):=n^?(^^^'^) (2.9) 

and ij is the j*^ component of the multi-index i. 

Now consider the tensor product of one-dimensional polynomial projection operators, each 
of degree n^: 

(2.10) 

This is an operator that projects elements of H*-*^-* onto a multivariate polynomial space whose 
truncation is given by the vector n — (ni, . . . , Ud) G Ng. This truncation is independent in each 
dimension, such that the range of the projector contains polynomials of degree up to Ui in each 
input x^^\ 

^i''^(/) = E---E(/*i)*i (2-11) 

il=0 id=0 

As in the one-dimensional case, the truncation induces error equal to the sum of the squares of 
the omitted coefficients, which we may similarly reduce to zero as ni — > oo,yi. The multivariate 
polynomial expansion also converges in an sense for any / £ 'HS'^\ 

3. One Dimensional and Full Tensor Approximations. The previous section intro- 
duced the two problems we explore throughout this paper: integration and projection onto 
polynomial spaces. Our next step is to develop useful (computable) approximation algorithms 
for these problems. To this end, we employ quadrature and pseudospectral approximation, re- 
spectively. Furthermore, we introduce criteria to describe when these approximations are precise. 

3.1. General setting. Given a one-dimensional linear operator we work with a con- 
vergent sequence of computable approximations, £;„ , such that 

^Oasm^cx) (3.1) 

in some appropriate norm. Taking the tensor product of these approximations provides an 
approximation to the full tensor operator: 

L^''^^C^:^^:^C^^\®---®Lt\. (3.2) 

The level of the approximation may be individually selected in each dimension, so the tensor ap- 
proximation is identified by a multi-index. As a consequence of the one-dimensional convergence, 
we have 



ll'C^'^^-r^^ll-^Oasm^oo. 



(3.3) 
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An important property of the approximation algorithms in which we are interested is whether 
they are precise for some inputs. We define the accurate set as the set of inputs for which the 
approximation is exact. 

Definition 3.1 (Accurate Sets). For an operator C and a corresponding approximation 
Cm, define the accurate set as A{Cm) '■— {/ : C{f) = C„i{f)}. 

It will also be useful (e.g., for quadrature) to consider when application of the operator to 
P is precise. 

Definition 3.2 (Half Accurate Sets). For an operator C and a corresponding approximation 
Cm, define the half accurate set A2{Cm) ■= {/ : C{p) = Cm{P)}- 

By linearity of the operator and its approximations, we can prove the following theorem 
relating the accurate sets of one-dimensional approximations and tensor approximations. 

Theorem 3.3. IJ a tensor approximation C^m is constructed from one- dimensional approx- 
imations Cm with known accurate sets, then 

A{C^^\) ® ■ ■ ■ <E> A{Clil) c A{C^^^) (3.4) 

Proof. It is sufficient to show that the approximation is precise for an arbitrary input /(x) = 
/(i)(a;(i))/(2)(a,(2)) . . . /(d)(a;(rf)) ^here /(^)(a;(')) G ^(-CmJ, because we may extend it to sums of 
products by linearity: 

^(d)(;(i)...^(.))^£(i)(^(i))^...^ 

= (s>---(s> c^'^\f^'^^) = c'^'^Xf). (3.5) 

This replacement is by linearity of tensor products and the definition of accurate sets, proving 
the statement. □ 

3.2. One-dimensional quadrature. Numerical quadrature approximates the action of an 
integral operator I^*^ with a weighted sum of point evaluations. For some family of quadrature 
rules, we write the "level m" quadrature rule, comprised of p^'^(m) : N -> N points, as 

p'*' (m) 

^^^H/) « s^;H/) := E ^ff^) (3-6) 

The choice of p('^(m) depends on the quadrature family, as some rules only exist for certain 
numbers of points. For others, we tailor p^^\m) depending on the desired properties, e.g., linear 
or exponential growth. 

Many quadrature families are precise if / is a polynomial of a particular degree or less, which 
allows us to specify a part of the accurate set. This definition leads to the following lemma. 
Lemma 3.4. For a one- dimensional quadrature rule with known polynomial accuracy a^'^\m), 

PaOM C ^(Q(^') and Pfloor(aW(™)/2) ^ ^2(2^^ 

It is intuitive and useful in later sections to draw the accurate set, as in the following example: 
Example 3.5. Figure 3.1 depicts the accurate set of a one- dimensional Gaussian quadrature 
rule with three points, Q^^\ which is accurate for polynomials up to fifth order. In this figure, the 
dots represent increasing polynomials orders of Xi. The solid box represents the extent of A{Q^^^), 
hence only polynomials inside or on the box are correctly integrated in general. The dashed box 
depicts the half accurate set, A2{Q''3^'^), drawn as half the geometric size. This skips the step of 
rounding down, but is visually more intuitive and does not change the result. 

Most quadrature families are constructed so that Qm -> l'*-* as m ^ oo for a wide variety 
of functions /. For this work, we rely on quadrature rules that exhibit polynomial accuracy of 
increasing order, which is sufficient to demonstrate convergence for functions in L^. 
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Accurate Set 

- - - Half Accurate Set 



012345678 
Order of polynomial 

Figure 3.1. The accurate set A{Q^) and half accurate set A2{Q^) of a three point Gaussian quadrature 
rule. Note that the vertical axis is meaningless here, as the figure is actually one dimensional. 

3.3. Multi-indices. Before continuing, we must make a short diversion to multi-indices, 
which provide helpful notation when dealing with tensor problems. A multi-index is a vector 
i e Nq. An important notion for multi-indices are neighborhoods. 

Definition 3.6 (Neighborhoods of multi- indices) . A forward neighborhood of a multi-index 
k is the multi-index set n/(k) := {k-|-ej : £ {1 . . . d}}, where Bi are the canonical unit vectors. 
The backward neighborhood of a multi-index k is the multi-index set nb(k) := {k — e, : Vi e 
{1...4,k-e, eNf?}. 

Smolyak algorithms rely on multi-index sets that are admissible. 

Definition 3.7 (Admissible multi-indices and multi-index sets). A multi-index k is admis- 
sible with respect to a multi-index set JC if nb{\s.) C IC. A multi-index set K is admissible if every 
Vi € IC is admissible with respect to IC. 

Example 3.8. The two dimensional index set, K c Ng, 

/C = {(0,0), (0,1), (1,0), (2,0)} (3.7) 

is admissible because the backward neighbors of all the multi-indices in the set IC are also in 
IC. For example, nf,((0, 0)) = {} because its backward neighbors are not actually in Nq. Also, 
nb((l,0)) = {(0,0)}, which is in the set. There are indices that are admissible with respect to 
IC, but not in IC: (3,0), (0,2), or (1,1). Any of these may be added to IC to produce another 
admissible index set. The multi-indices (3, 1), (4, 0), or (4, 4) are inadmissible with respect to IC, 
so adding them to IC would produce an inadmissible m/alti-index set. 

Two common admissible multi-index sets with simple geometric structure are total order 
multi-index sets and full tensor multi-index sets. We often encounter total order sets in the 
sparse grids literature and full tensor sets when dealing with tensor grids of points. The total 
order multi-index set /C^ comprises those multi-indices that lie within a d-dimensional simplex 
of side length n: 

/C^={kGNM|k||i<n} (3.8) 

The full tensor multi-index set /Cn is the complete grid of indices bounded term-wise by a multi- 
index n; 

/CJ; := {k e : Vi e {1 . . . d}, h < m) (3.9) 

3.4. Tensor product quadrature. Using the linearity of integral operators, we can con- 
struct the full tensor quadrature approximation to a multi-dimensional integral I^**) by taking 
the tensor product of quadrature operators: 
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2 4 6 

Order of polynomial 



Figure 3.2. The quadrature rule accurate set A{Q)J^,) and half accurate set A2(QtJ.,^) for a two dime 
sional rule constructed from three point Gaussian quadrature rules. 



where 



7 — IC^ 



(3.11) 
(3.12) 

(3.13) 



This quadrature rule is completely defined by the choice of one-dimensional quadrature families, 
and the nodes and weights are just combinations of the constituent one-dimensional rules. As a 
consequence of Theorem 3.3, we can describe the accurate set of the tensor algorithm. 

Corollary 3.9. The accurate set of a tensor quadrature rule constructed from constituent 



approximations with known accurate sets includes yl(Qmi) 



'^(Qi'^) c^(Q^)). 



(2) 

Example 3.10. Let 8(3 3) be the two-dimensional full tensor quadrature rule constructed 



from three point Gaussian quadrature rules. Then 
1(2 



3,(2) 



C -4(Q^3 3j) and iL 2 ^'^2 

^2(2(33))- These accurate sets are depicted in Figure 3.2. 

While Corollary 3.9 specifies the natural extension of the accurate set to the tensor case, the 
accurate set may also include other functions that have more structure: 

Lemma 3.11. Let f be separable, without loss of generality, with respect to its first input, 
i.e., /(x) = f{xi)f'{x2,...Xd). Further, let X^{f') = and let /'(xi) G yl(Q^l). Then 
Qia if) — T^'^^f) — for any /". Hence, / G yl(Q^^). This is a useful property arising from 
the linearity of integrals, which we use in our analysis of sparse problems. 

By the same argument as in the one-dimensional case, for functions / G we have 

Q^H/)->l('^H/)- 

3.5. Pseudospectral approximation in one dimension. Pseudospectral approximation 
provides a practical non-intrusive algorithm by approximating a truncated polynomial chaos 
expansion with quadrature operators. Define the pseudospectral approximation in one dimension 
as 



j(2) 



c 



q(''(m) 
q^^^ (m) 



(3.14) 
(3.15) 
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where q^^\m) is the polynomial truncation, to be specified shortly [4, 16]. Pseudospectral ap- 
proximations are constructed around a level m quadrature rule, and may include as many terms 
in the sum as possible while maintaining accuracy. Assuming that f G L'^, we can compute 
the error between the pseudospectral approximation and an exact projection onto the same 
polynomial space: 



(j^*' (m) g*-*' (m) 

i=o k=0 



= E (/. - 

3=0 



(3.16) 



This quantity is the aliasing error [4, 16]. The error is non-zero because quadrature in general 
only approximates integrals, hence each fi is an approximation of fi. The pseudospectral operator 
also incurs truncation error, as before, which is orthogonal to the aliasing error. We can expand 
each approximate coefficient 



= E/.Q«(V'?V«) 



k=0 



oo 



A;=0 



/=q(i)(m) + l 



(3.17) 
(3.18) 

(3.19) 



The first step substitutes in the polynomial expansion of /, which we assume is convergent, 
and rearranges using linearity. The second step partitions the sum around the truncation of 
the pseudospectral expansion. Although V'i*'') = Sjk, we cannot assume in general that 



Q 



Now substitute (3.19) back into the aliasing error expression: 



'(m) 



'(m) 



q'-''\m) oo 

E E E E 

j=0 j=0 \ fe=0 /=g(')(m) + l 

(3.20) 

This form reveals the intimate link between the accuracy of pseudospectral approximations 
and the polynomial accuracy of quadrature rules. All aliasing is attributed to the inability of the 
quadrature rule to determine the orthonormality of the basis polynomials, causing one coefficient 
to corrupt another. The contribution of the first two parenthetical terms on the right of (3.20) is 
called internal aliasing while the third term is called external aliasing. Internal aliasing is due to 



)(») 



inaccuracies in Q(V'|*''^fc*'') when both ip'",''' and ■0^''' are included in the expansion, while external 



(i) 



aliasing occurs when only one of these polynomials is included in the expansion 
practical quadrature rules — and for all those used in this work 



whenever (V'^^^V'!*'') 



For many 
-4(2) 



and hence the discrete inner product is not zero, Q{ipj^''')pk^) is [34]. This provides an 

important or der-of- magnitude estimate of the error. 

Both types of aliasing errors are driven to zero by sufficiently powerful quadrature, but we 
are left with a practical trade-off. The obvious option is to select the smallest quadrature level 
that sets the internal aliasing to zero, i.e., that allows the pseudospectral operator to precisely 
recover functions that actually lie in the space spanned by polynomials included in the truncated 
expansion. We refer to this space as the range of the pseudospectral projection operator. 

Furthermore, recall that because ^ ff < oo, the high-order coefficients must in general 
decay. Aliasing error is proportional to coefficient magnitude because ||Q('0i'''V'i*'')ll is 
Therefore, we expect that for a pseudospectral method to accurately approximate a particular 
/, the possible errors resulting from internal aliasing are much larger than those from external 
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aliasing. Indeed, we see in practice that any internal aliasing can lead to hopeless corruption of 
the approximation. 

Analysis of (3.20) yields a theorem defining when the pseudospectral projection operator has 
no internal aliasing and thus is precise on its range space: 

Theorem 3.12. The internal aliasing of is zero if^j, k < q^'-^m), {ipfip'^!^) G ^(Qm ), 
or equivalently, Vj < ^(^^(to), G ^2(Qm )• Thus Pj^,)(„) C A{S^!^). 

The most efficient choice is to include as many basis terms as possible without introducing 
internal aliasing error. 

Lemma 3.13. The conditions of Theorem 3.12 are satisfied if q^'^\m) — floor(a^*-'(m)/2) 
when Qm has polynomial accuracy a*^*-'(rn). 

Example 3.14. Revisit Figure 3.1, which shows accurate sets for a three-point Gaussian 
quadrature rule. The dashed box, -42(23), shows all the terms that can be included in the 
pseudospectral approximation constructed from this quadrature rule without inducing internal 
aliasing — that is, terms which are precisely recovered by the operator. In general, an m-point 
Gaussian quadrature rule allows one to compute polynomial approximations up to order m — 1. 

Given that Theorem 3.12 is satisfied, we wish to show that the pseudospectral approximation 
converges to the true function, where the magnitude of the error is as follows: 



00 \ 00 



/-5«(/)|[= E E hQi^f^^')] + E f?- (3-21) 

j=0 \fe=<j(')(m) + l / ;=g(»)(m) + l 

The two terms on right hand side comprise the external aliasing and the truncation error, respec- 
tively. We already know that the truncation error goes to zero as q^^\m) — > 00. The external 
aliasing also vanishes for functions / e as the truncated portion of / likewise decreases [34]. 
In the case of Gaussian quadrature rules, a link to interpolation provides precise rates for the 
convergence of the pseudospectral operator based on the regularity of / [4] . 

3.6. Full tensor pseudospectral approximation in d dimensions. As with quadrature 
algorithms, pseudospectral approximation in one dimension is directly extensible to full tensor 
problems by taking the tensor product of the one-dimensional approximations, 

S^^:^S^l^---^S^l^ E Q^'(/^k)^k(x) (3.22) 

where q(m) — ((/'-^'(toi), . . . ,q^'^\mj)). This expansion contains all basis functions up to and 
including ^'q(m)- As in the one-dimensional case, error in the approximation can be separated 
into truncation, internal aliasing, and external aliasing. We may again apply Theorem 3.3 to 
derive the accuracy of the tensor product algorithm. 

Corollary 3.15. The accurate set of a tensor product pseudospectral approximation oper- 
ator constructed from constituent approximations with known accurate sets includes A{Sm\) ® 

We conclude that the full tensor approximation has no internal aliasing. 

Example 3.16. Figure 3.2 shows the accurate and half-accurate sets of Q^gg-j, indicating 
which terms may be included in the pseudospectral expansion without introducing internal aliasing 
error. 

The same analysis as in the one-dimensional case allows us to conclude that the d-dimensional 
pseudospectral approximation converges to the true function. 

4. Smolyak Algorithms. Thus far, we have developed polynomial approximations of 
multivariate functions by taking tensor products of one-dimensional pseudospectral operators. 
Smolyak algorithms avoid the exponential cost of full tensor products by assuming that the input 
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dimensions are not fully coupled, and then using a telescoping sum to blend different full tensor 
approximations . 

Example 4.1. Assume that f{x,y) = + + x^y. To construct a polynomial expansion 
with both the and y"^ terms, a full tensor pseudo spectral algorithm would estimate all the 
polynomial terms up to x'^y'^ , because tensor algorithms fully couple the dimensions. This mixed 
term is costly, requiring, for example, an 8x8 point grid for Gaussian quadratures. The individual 
terms can be had much more cheaply, using 8x1, 1x8, and 4x2 grids, respectively. Smolyak 
algorithms help realize such savings in practice. 

4.1. General Smolyak algorithms. As in Section 3, assume that we have for every di- 
mension i = 1 ... d a convergent sequence >C^*^ of approximations to a one-dimensional linear 
operator such that \\C^^^ - d{^\\ as — > oo. Let C denote the collection of these 

sequences over all the dimensions. For any i, we may write the exact or "true" operator as the 
telescoping series 



(4.1) 

fe=0 



Define the difference operators 



A« := Cf = 0, (4.2) 
AW:=/:W-/:«i. (4.3) 

These allow us to rewrite the telescoping sum as 



CS^^Y.^- (4-4) 

fc=0 



Now we may write the tensor product of the exact operators as the tensor product of the tele- 
scoping sums, and interchange the product and sum: 



oo 



oo 

= ^A(.>...«Ai^j (4.6) 

k=0 

Smolyak's idea is to approximate the tensor product operator with truncations of this sum [31]: 

^(/C,d,/:):=^Ai>...®A(J. (4.7) 

kG/C 

We refer to the multi-index set K- as the Smolyak multi-index set, and it must be admissible for 
the sum to telescope correctly. Smolyak specifically suggested truncating with a total order multi- 
index set, which is the most widely studied choice. However, we can compute the approximation 
with any admissible multi-index set. Although the expression above is especially clean, it is 
not the most useful form for computation. We can reorganize the terms of (4.7) to construct a 
weighted sum of the tensor operators: 

A{JC, d,C)^Yl ^k^i? • • • C^k} ' (4-8) 

kGK 



where Ck are integer Smolyak coefficients computed from the combinatorics of the difference 
formulation. One can compute the coefficients through a simple iteration over the index set and 
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use (4.7) to determine which full tensor rules are incremented or decremented. In general, these 
coefficients are non-zero near the leading surface of the Smolyak multi-index set, reflecting the 
mixing of the most accurate constituent full tensor approximations. 

If each sequence of one-dimensional operators converges, then the Smolyak approximation 
converges to the tensor product of exact operators as /C H> Nq. For the isotropic simplex index 
set, some precise rates of convergence are known with respect to the side length of the simplex 
[36, 37, 35, 28, 30]. Although general admissible Smolyak multi-index sets are difficult to study 
theoretically, they allow detailed customization to the anisotropy of a particular function. 

4.2. Accuracy of Smolyak algorithms. In the one-dimensional and full tensor settings, 
we have characterized approximation algorithms through their accurate sets — those inputs for 
which the algorithm is precise. This section shows that if the constituent one-dimensional ap- 
proximations have nested accurate sets, Smolyak algorithms are the ideal blending of different 
full tensor approximations from the perspective of accurate sets; that is, the accurate set of the 
Smolyak algorithm contains the union of the accurate sets of the component full tensor approxi- 
mations. This result will facilitate subsequent analysis of sparse quadrature and pseudospectral 
approximation algorithms. 

Theorem 4.2. Let A(IC,d,C) be a Smolyak algorithm composed of linear operators with 
nested accurate sets. That is, let K, he admissible and m < m! imply A{Cm) ^ A{c''^,) for 
i — 1 . . . d. Then the accurate set of A(IC, d, C) contains 

A (A(/C, d, £)) D y ^ (r^f • • • ® £[^J) (4.9) 

D y ^(£W)®...®^(4'j). (4.10) 

ke/c 

Our proof closely follows the framework provided by Novak and Ritter, but includes a gen- 
eralization to arbitrary Smolyak multi-index sets [24, 25, 2]. We begin by introducing notation 
to incrementally build a multi-index set dimension by dimension. 

Definition 4.3. For a multi-index set K, of dimension d, let the restriction of the multi- 
indices to the first i dimensions be /C'*' := {ki:i = {ki, . . . ,ki) : k G JC\. Furthermore, define 
subsets of K, based on the i*^ element of the multi-indices, /cj'^ := {ki:i : k £ /C*^*' and ki^i — j}. 
These sets are nested, /cj'^ 3 ICj^^i, because K, is admissible. Also let kf^'^^ denote the maximum 
value of the i*^ component of the multi-indices in the set K.. 

Using this notation, one can construct K. inductively, 

/C(i' = {l,...,/cr''} 
/CW = y icl'-''>^j, 1^2. ..d 

Proof of Theorem 4.2. It is sufficient to prove that the Smolyak operator is precise for 
an arbitrary / with tensor structure, f — fi ® ■ ■ ■ ® fd- Suppose there exists a k* such that 
/ e A{C''^}). We will show that if JC is an admissible multi-index set containing k*, then 
A{lC,d,C) is precise on /. We do so by induction on the dimension i of the Smolyak operator 
and the function. 

First, consider the i = 1 case. A{JC^^^ , 1, -C) = /:iiL , where > kl . Hence A{A{JC^^^ , 1, -C)) = 

^(rilL). 

For the induction step, we construct the (i + l)-dimensional Smolyak operator in terms of 
the i-dimensional operator: 

^(/C(''+i),z+ 1,/:) = A(/Cf ,»,£) ® (4^+^) - C!^+^). (4.13) 



(4.11) 
(4.12) 
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This sum is over increasing levels of accuracy in the i+1 dimension. We know the level required 
for the approximate operator to be precise in this dimension; this may be expressed as 

4+'\f.+i) = 4^,'\f^+i) = C^'+'H^^+l) when J - 1 > k*^,. (4.14) 

Therefore the sum (4.13) can be truncated at the k*_^-^ term, as the differences of higher terms 
are zero when applied to /: 

A{JC^^+'\i +1,C) = J2 C) ® (4'+'^ - (4.15) 

Naturally, k^.^ e idfl . By nestedness, kj.^ is also contained in /cj*'' for j < k*^-^^. The induction 
hypothesis then guarantees 

/i ® • • • ® /, e A{A{lcf\t,C)), Vj < fc*+i. (4.16) 
Applying the (i + l)-dimensional Smolyak operator to the truncated version of / yields 

= J2 Aijcf.i, C)Ui ®---®f.)® (4+'^ - 4_+/^)(/.+i). (4.17) 

3 = 1 

Since each of the i-dimensional Smolyak algorithms is precise, by the induction hypothesis, we 
replace them with the true operators and rearrange by linearity to obtain 

A(/C('+i),* + l,/:)(/i /.+i) - £(■)(/! ®---®f.)® 5^(4'+'^ - 4^_\^))(/,+i)(4.18) 

= £(')(/! ®---®fi)® 4^|)(/,+i). (4.19) 

The approximation in the i + 1 dimension is exactly of the level needed to be precise on the 
{i + l)*'^ component of /. Then (4.19) becomes 

£«(/! C^'+^HUi) = >C('+^H/i ® • • • ® /.+i) (4.20) 

Thus the Smolyak operator is precise for /, and the claim is proven. □ 

4.3. Smolyak quadrature. We recall the most familiar use of Smolyak algorithms, sparse 
quadrature. Following the formalism of Section 4.1, consider a family of one-dimensional quadra- 
ture rules ^Q^*''^ in each dimension i = I . . .d; denote these rules by Q. We define Qq*'' :— 
as required for the Smolyak construction. A Smolyak quadrature algorithm in d dimensions is a 
weighted sum of full tensor quadrature rules: 

^(/C,d,Q)-^CkQi''^ (4.21) 

The set of functions that are integrated correctly by a Smolyak quadrature algorithm is described 
as a corollary of Theorem 4.2. 

Corollary 4.4. For a sparse quadrature rule satisfying the hypotheses of Theorem 4-2 (that 
is, composed of one- dimensional quadrature rules that have nested accurate sets), we have 

A{A{IC,d,Q))D [J A{Q^^^) (4.22) 
ke/c 
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Figure 4.1. The accurate set diagram for two Smolyak quadrature rules, and the corresponding basis for a 
Smolyak pseudospectral estimate. ^(Q) is shown in a solid line, ^2(Q) is the dashed line. 
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Additionally, other separable functions may be precisely integrated. 

Theorem 4.5. /// is separable such that f — fi{xi)f'{x2, • • ■ , a^d) and X^^\fi) — 0, and if 
for every full tensor rule k G /C we have fi £ A{Q\}^), then f G A{A{IC, d, Q)). 

Proof. If every full tensor rule in the Smolyak algorithm correctly computes that /i inte- 
grates to zero, then their weighted sum is necessarily zero. □ 

We conclude by illustrating accurate sets for some sparse quadrature rules. 

Example 4.6. Figure 4-.1 shows the accurate set for two instances of Smolyak quadrature, 
constructed from one dimensional Gaussian rules with different growth rates p(^'(m). The accu- 
rate sets look like steps, as they are the superposition of rectangular full tensor accurate sets. 

4.4. Smolyak pseudospectral approximation. Applying Smolyak's algorithm to pseu- 
dospectral approximation operators yields a sparse algorithm that converges under similar con- 
ditions as the one-dimensional operators it is constructed from. This algorithm is written 

A{IC,d,S)^J2^^^k^- (4-23) 

Remark 4.7. The Smolyak algorithm is a sum of different full tensor pseudospectral approx- 
imations, where each approximation is built around the polynomial accuracy of a single full tensor 
quadrature rule. It is not naturally expressed as a set of formulas for the polynomial coefficients, 
because different approximations include different polynomials. 

Lemma 4.8. When the one- dimensional rules are constructed to satisfy Lemma 3.13, the 
term ^'j is included in the Smolyak approximation if and only i/ 3k e /C : ^'j G ^2(Qk*'')- Here, 
is a full tensor quadrature rule corresponding to the quadrature rule used by the full tensor 
pseudospectral approximation sj^^ . 

Example 4.9. Returning to Figure 4-1, the half accurate set again depicts the polynomial 
terms included in the corresponding Smolyak pseudospectral approximation. 

In principle, we could stop our analysis here, as the convergence of the sparse algorithm 
is assured by the properties of Smolyak's algorithm. However, for the purposes of comparing 
to direct quadrature, we can provide two theorems that characterize the internal aliasing and 
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external aliasing in more detail. Specific examples of internal and external behavior are given in 
the next section. The internal aliasing result follows as a corollary of Theorem 4.2. 

Corollary 4.10. // the constituent one-dimensional pseudo spectral rules are consistent 
with Theorem 3.12, then the resulting Smolyak pseudospectral algorithm has no internal aliasing 
error. In other words, the Smolyak pseudospectral approximation is precise on its range. 

Theorem 4.11. Let ^'j he a polynomial term included in the expansion provided by the 
Smolyak algorithm A(IC,d,S), and let ^>y be a polynomial term not included in the expansion. 
There is no external aliasing of "^y onto ^'j if either (a) for any dimension i, j[ < ji or (b) 
3k G /C such that is included in the expansion provided by S^^^ and ^'j'^'j G A{Q^^^), where 
Q corresponds to the quadrature rules of S . 

Proof. If condition (a) is satisfied, then ^'j and "^y are orthogonal in dimension i, and hence 
that inner product is zero. Every quadrature rule that computes the coefficient of is accurate 
for polynomials of at least order 2j. Since j[ + ji < 2ji, every rule that computes the coefficient 
can determine the orthogonality, and therefore there is no aliasing. If condition (b) is satisfied, 
then the result follows from the cancellations the Smolyak algorithm exploits, as seen in the proof 
of Theorem 4.2. □ 

The properties stated in these theorems are extremely useful. First, any Smolyak pseu- 
dospectral algorithm, regardless of the admissible Smolyak multi-index set used, has no internal 
aliasing. Second, there is external afiasing, as we expect, but the algorithm uses basis orthogonal- 
ity to limit which external coefficients can alias onto an included coefficient. Hence, the Smolyak 
construction is a useful approximation, in that we can tailor it to perform the desired amount 
of work and achieve reliable approximations of the selected coefficients. Computing an accurate 
approximation of the function only requires including enough terms so that the truncation and 
external aliasing errors are small. 

5. Comparing Direct Quadrature to Smolyak Pseudospectral Approximation. 

The current UQ literature often suggests a direct quadrature approach for constructing polynomial 
chaos expansions [40, 20, 7, 17]. In this section, we describe this approach and show that, in 
comparison to a true Smolyak algorithm, it is inaccurate or inefficient in almost all cases. 

5.1. Direct quadrature polynomial expansions. At first glance, direct quadrature is 
quite simple. First, choose a multi-index set J7 to define a truncated polynomial expansion: 



The index set is typically admissible, but need not be. Second, select any d-dimensional 
quadrature rule Q^'^\ and estimate every coefficient as: 



Unlike the Smolyak approach, we are left to choose and Q^**' independently, giving the appear- 
ance of ffexibility. We are interested in selecting Q and J' identically to the Smolyak approach, 
as in Lemma 4.8. 

5.2. Internal aliasing in direct quadrature. As before, we define internal aliasing as 
the error incurred when one coefficient in the expansion aliases onto another because the included 
basis functions are not numerically orthogonal. 

Theorem 5.1. For a multi-index set J and a quadrature rule Q}''^\ the corresponding direct 
quadrature polynomial expansion has no internal aliasing if and only if for every j,y G J , ^j^j' 
is correctly integrable, that is, lies in A{Q}''^^\ 

This theorem follows directly from the definition of accurate sets. We can immediately 
conclude that for any basis set J , there is some quadrature rule sufficiently powerful to avoid 
internal aliasing error. In practice, however, this rule may not be a desirable one. 

Example 5.2. Assume that for some function with two inputs, we wish to include the 
polynomial basis terms (a,0) and (0, &). By Theorem 5.1, the product of these two terms must 




(5.1) 



(5.2) 
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be in the accurate set; hence, the quadrature must include at least a full tensor rule of accuracy 
(a, 6). Although we have not asked for any coupling, direct quadrature must assume full coupling 
of the problem in order to avoid internal aliasing. 

Thus we reach the surprising conclusion that direct quadrature is essentially incompatible 
with actually using sparse grid evaluations. For most sparse quadrature rules, we cannot include 
as many terms as the Smolyak pseudospectral approach without incurring internal aliasing, be- 
cause the quadrature is not powerful enough in the mixed dimensions. 

Example 5.3. Lei X be the two-dimensional domain [—1, 1]^. Select a uniform weight func- 
tion, which corresponds to Legendre polynomials and Gauss-Legendre quadrature. Let f{x, y) = 
'4>o{x)'4>4{y) ■ Use an exponential growth rule for the quadrature, such that p^'^\m) — 2™^^. 5*6- 
lect a sparse quadrature rule based on a total order multi-index set Qi-t . Figure 5.1 shows the 

accurate set of this Smolyak quadrature rule (solid line) along with its half-accurate set (dashed 
line), which encompasses all the terms in the direct quadrature PCE. 

Consider the j = (8, 0) polynomial, which is in the half-accurate set. The product of the 
(0,4) and (8,0) polynomial terms is (8,4), which is not within the accurate set of the sparse rule. 
Hence, (0,4) aliases onto (8,0) because this quadrature rule has limited accuracy in the mixed 
dimensions. 

Using both the Smolyak pseudospectral and direct quadrature methods, we numerically com- 
puted the polynomial expansion for this example. The resulting coefficients are shown in Figure 
5.2. Even though the two methods use the same information and project f onto the same ba- 
sis, the Smolyak result has no internal aliasing while direct quadrature shows significant internal 
aliasing. Although both methods correctly compute the (0, 4) coefficient, direct quadrature shows 
aliasing on (8,0) as predicted, and also on (10,0), (12,0), and (14,0). This polynomial approx- 
imation is so drastically corrupted as to be completely useless. Alternating terms are computed 
correctly because of the parity of the functions. 

Figure 5.1. The accurate set and polynomials included in the direct quadrature construction from Example 5.3. 
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This example depicts a dangerous type of internal aliasing: a low order term, which we expect 
to be large, corrupts many other terms, which we would ordinarily expect to be much smaller. 

We can understand this phenomena geometrically, in that whenever the accurate set A is 
concave by more than unit steps, then the product of some pair of polynomials in A2 (drawn as 
the sum of the points), does not lie in A. Figure 5.1 shows that A is quite concave, leading to the 
observed errors. Forbidding all sparse quadrature rules with concave accurate sets — which imply 
very weak coupling between inputs — is quite restrictive! This generally excludes every nested 
quadrature rule, or any super-linear growth in the number of points per level. 

Remark 5.4. Sparse quadrature rules constructed from Gaussian quadrature rules with 
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Figure 5.2. Numerical results for Example 5.3; each color square indicates the log of the coefficient magni- 
tude for the basis function at that position. The circle identifies the correct non-zero coefficient. 




Order of polynomial Order of polynomial 

(a) Smolyak Pseudospectral (b) Direct Quadrature 

p'^'^\m) = m, truncated according to an isotropic total-order multi-index set, do not have internal 
aliasing. Although the accurate set has a stair-step shape of size two, this shape is "convex" 
enough. 

The results and analysis in this section underpin one of key conclusions of this work: direct 
quadrature strategies fail to accurately compute polynomial expansions for most sparse quadra- 
ture rules, if one truncates the polynomial basis similarly to the Smolyak approach. It is possible 
to select a sufficiently small polynomial basis to avoid internal aliasing, but this approach requires 
drastic conservatism that could easily be avoided with a Smolyak pseudospectral approximation. 

5.3. External aliasing. The difference in external aliasing between direct quadrature and 
Smolyak pseudospectral approximation is much less severe. Both methods exhibit external alias- 
ing from terms far outside the estimated set. And since they are constructed from similar con- 
stituent one-dimensional quadrature rules, the aliasing is of similar magnitude when it occurs. 
The fundamental theorem is quite similar to Theorem 4.11. 

Theorem 5.5. For a multi-index set J and a quadrature rule Q^'^\ the corresponding 
direct quadrature PCE approximation has no external aliasing between an included term ^Pj and 
a non-included term if^y^j G A{Q^'^''). 

Yet the two methods perform differently because of their behavior on separable functions. 
The Smolyak pseudospectral algorithm gains a condition that prevents external aliasing as a 
result of Theorem 4.5 and thus it has strictly less external aliasing in general. 

Example 5.6. // we repeat Example 5.3 hut choose f to he a polynomial outside the approxi- 
mation hasis, f = '06(2;)V'6(y); we ohtain the results in Figure 5.3. Now every non-zero coefficient 
is the result of external aliasing. The Smolyak pseudospectral approach incorrectly constructs sev- 
eral coefficients, hut aliasing occurs only for terms of lower order, hetween (0,0) and (6,6). Some 
terms are integrated correctly hecause of parity or hecause (6,6) is sufficiently close to the approx- 
imation hasis for the terms to satisfy the criterion of Theorem -<^.ii, ^'/^'j G ^(Qk*"*) for some 
k e /C. 

Direct quadrature makes the same mistakes, hut is not limited to errors helow (6, 6) hecause 
Theorem 5.5 has different conditions; hence the aliasing is more severe. 

This example is representative of the general case. Direct quadrature always incurs at least 
as much external aliasing as the Smolyak approach, and the methods become equivalent if the 
external term causing aliasing is of very high order. Although both methods exhibit external 
aliasing onto coefficients of the approximating basis (the "internal coefficients"), in a realistic 
problem (with suitable decay and truncation) the magnitude of the external coefficients and 
the corresponding corruption should be much smaller than the internal coefficients themselves. 
Hence, the error introduced in the approximation is manageable. 
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Figure 5.3. Numerical results for Example 5.6; each color square indicates the log of the coefficient mag- 
nitude for the basis function at its position. The circle indicates the correct non-zero coefficient. The Smolyak 
pseudospectral approach has fewer terms corrupted by external aliasing in this case. 
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5.4. Summary of comparison. Compared to the Smolyak pseudospectral approach, di- 
rect quadrature yields larger internal and external aliasing errors. Because of these aliasing 
errors, direct quadrature is essentially unable to make efficient use of most sparse quadrature 
rules. The Smolyak pseudospectral approach, on the other hand, is guaranteed never to have 
internal aliasing if the one-dimensional pseudospectral operators are chosen according to simple 
guidelines. We therefore recommend against using direct quadrature. 

6. Adaptive Polynomial Expansions. When constructing a polynomial approximation 
of a black-box computational model, there are two essential questions: first, which basis terms 
should be included in the expansion; and second, what are the coefficients of those basis terms? 
The Smolyak construction allows detailed control over the truncation of the polynomial expansion 
and the work required to compute it. Since we typically do not have a priori information about 
the functional relationship generated by a black-box model, we develop an adaptive approach to 
tailor the Smolyak approximation to this function, following the dimension-adaptive quadrature 
approach of Gerstner & Griebel [13] . 

The Smolyak algorithm is well suited to an adaptive approach. The telescoping sum converges 
as the index set grows, so we can simply add more terms to improve the approximation until we 
are satisfied. We separate our adaptive approach into two components: local improvements to 
the Smolyak multi-index set and a global stopping criteria. 

6.1. Local adaptivity. Local adaptivity is responsible for identifying multi-indices to add 
to a Smolyak multi-index set in order to improve its accuracy. A standard approach is simply to 
grow an isotropic simplex of side length n. Gerstner & Griebel instead suggest a series of greedy 
refinements that customize the Smolyak algorithm to a particular problem [13]. 

The local refinement used in [13] is to select a multi-index k S /C and to add the forward 
neighbors of k that are admissible. The multi-index k is selected via an error indicator e(k). 
We follow [13] and assume that whenever k has a large impact on the result of the algorithm, it 
represents a "direction" that likely needs further refinement. 

Let k be a multi-index such that /C' := /C U k, where K. and JC' are admissible multi-index 
sets. The triangle inequality (for some appropriate norm) bounds the change in the Smolyak 
approximation produced by adding k to /C, yielding a useful error indicator: 

\\A{K\d,C)^A{K,d,C)\\ < ||Ai^®-..®Atl| -: e(k) (6.1) 

Conveniently, this error indicator does not change as K. evolves, so we need only compute it once. 
At each adaptation step, we find the k that maximizes e(k) and that has at least one admissible 
forward neighbor. Then we add those forward neighbors. 
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6.2. Global adaptivity. Now that we have a strategy to locally improve a multi-index set, 
it is useful to have a global estimate of the error of the approximation, eg. We cannot expect to 
compute the exact error, but even an indicator which is correlated with the error is useful. We 
follow Gerstner & Griebel's choice of global error indicator 



where the sum is taken over all the multi-indices that are eligible for local adaptation at any 
particular step (i.e., that have admissible forward neighbors). 

6.3. Comments. Thus far we have presented the adaptation strategy without reference to 
the problem of polynomial approximation. In this specific context, we use the i^(X, w) norm in 
(6.1) because it corresponds to the convergence properties of pseudospectral approximation. 

We conclude this section with a few observations about the adaptive approach. First, it 
is a heuristic method and we can only recommend it based on the success of adaptive sparse 
quadrature and experimental validation (see Section 9). Second, the greedy approach assumes 
that following the largest local errors leads to the largest global sources of error. Our analysis of 
external aliasing suggests that in the case of pseudospectral approximation, significant missing 
terms alias onto some of the included lower-order coefficients, giving the algorithm a useful 
indication of the direction to refine. Finally, Gerstner & Griebel suggest augmenting the local 
refinement algorithm to add terms that are computationally inexpensive, even if they do not 
appear immediately useful. In general, balancing the error indicator (6.1) with an indicator 
of the computational cost associated with a multi-index k might improve efficiency for certain 
problems, or avoid pathological cases where e is otherwise locally zero, but we do not explore 
this extension here. 

7. Selection of Quadrature Rules. Thus far we have sidestepped practical questions 
about which quadrature rules exist or are most efficient. Our analysis has relied only on polyno- 
mial accuracy of quadrature rules; all quadrature rules with a given polynomial accuracy allow 
the same truncation of a pseudospectral approximation. In practice, however, we care about 
the cumulative cost of the adaptive algorithm, which must step through successive levels of 
refinement. 

Integration over a bounded interval with uniform weighting offers the widest variety of 
quadrature choices, and thus allows a thorough discussion. Table 7.1 summarizes the costs 
of several common quadrature schemes. First, we see that linear-growth Gaussian quadrature is 
asymptotically much less efficient than exponential-growth in reaching any particular accuracy. 
However, for rules with fewer than about ten points, this difference is not yet significant. Second, 
Clenshaw-Curtis shows efficiency equivalent to exponential-growth Gaussian: both use n points 
to reach nth order polynomial accuracy [5] . However, their performance with respect to external 
aliasing differs: Clenshaw-Curtis slowly loses accuracy if the integrand is of order greater than n, 
while Gaussian quadrature gives 0(1) error even on (n-l- l)-order functions [34]. This may make 
Clenshaw-Curtis Smolyak pseudospectral estimates more efficient. Finally, we consider Gauss- 
Patterson quadrature, which is nested and has significantly higher polynomial accuracy — for a 
given cumulative cost — than the other types [26] . Computing the quadrature points and weights 
in finite precision (even extended-precision) arithmetic has practically limited Gauss-Patterson 
rules to 255 points, but we recommend them whenever this is sufficient. 

For most other weights and intervals, there are fewer choices that provide precise polynomial 
accuracy, so exponential-growth Gaussian quadrature is our default choice. In the specific case of 
Gaussian weight, Genz has provided a family of Kronod extensions, similar to Gauss-Patterson 
quadrature, which may be a useful option [12]. 

Remark 7.1. If a linear growth rule is chosen and the domain is symmetric, we suggest 
that each new level include at least two points, so that the corresponding basis grows by at least 
one even and one odd basis function. This removes the possibility for unexpected effects on the 
adaptive strategy if the target function is actually even or odd. 
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Table 7.1 

The cost of four quadrature strategies as their order increases: linear growth Gauss-Legendre quadrature (Lin. 
G), exponential growth Gauss-Legendre quadrature (Exp. G), Clenshaw-Gurtis quadrature (C-C), and Gauss- 
Patterson quadrature (G-P). We list the number of points used to compute the given rule (p), the polynomial 
accuracy (a), and the total number of points used so far (t). For nested rule, (p) = (t), so the total column is 
omitted. 



8. Connections to Interpolation. It is well known that in one dimension, pseudospec- 
tral approximation based on Gaussian quadrature is equivalent to Lagrange interpolation on 
the same quadrature points [4, 16, 3]. This connection to interpolation is used to prove the 
convergence of pseudospectral operators. Similarly, the tensor product of pseudospectral approx- 
imations based on Gaussian quadrature is equivalent to Lagrange interpolation on the full tensor 
grid. Unfortunately, the Smolyak pseudospectral algorithm we have described does not actually 
produce an interpolant, because none of the quadrature rules is both nested (Clenshaw-Curtis, 
Gauss-Patterson) and yields an interpolant (Gaussian quadrature). 

There is an extension to Lagrange interpolation called sparse interpolation, which uses in- 
formation from a sparse grid and has strong known convergence properties [23, 2]. Constantine 
et al. [6] show that the Smolyak pseudospectral method using Gaussian quadrature is equivalent 
to sparse interpolation. It is important to note that sparse interpolants are not necessarily inter- 
polants at every point on the sparse grid; instead, they are Smolyak combinations of full tensor 
interpolants. If the constituent full tensor rules are interpolants on nested point sets, then the 
Smolyak algorithm produces an interpolant on the corresponding sparse grid [29] . 

The connection of Gaussian quadrature-based pseudospectral operators to interpolation is 
useful for analysis, but it is an open question whether these rules are actually more efficient in 
practice, as compared to other options — specifically Clenshaw-Curtis or Gauss-Patterson. 

9. Numerical Experiments. Our numerical experiments focus on evaluating the perfor- 
mance of different quadrature rules and of the adaptive Smolyak approximation strategy. Aside 
from the numerical examples of Section 5, we do not investigate the performance of direct quadra- 
ture any further. Given our theoretical analysis of aliasing errors and the numerical demonstra- 
tions of Constantine et al., one can conclude without further demonstration that destructive 
internal aliasing indeed appears in practice. 

We begin by evaluating mean-square convergence on the Genz test functions. Then we 
approximate a chemical kinetic system, to illustrate the efficiency and accuracy of adaptive 
Smolyak pseudospectral methods. 

9.1. Basic convergence: Genz functions. The Genz family comprises six parameterized 
functions, defined on [— 1, 1]"^ — > M [10, 11]. They are commonly used to investigate the accuracy 
of quadrature rules and interpolation schemes [2, 19]. The purpose of this example is to show 
that different Smolyak pseudospectral strategies behave roughly as expected, as evidenced by 
decreasing approximation errors as more function evaluations are employed. 

Our first test uses five isotropic and non-adaptive pseudospectral approximation strategies. 
The first strategy is the isotropic full tensor pseudospectral algorithm, based on Gauss-Legendre 
quadrature, where the order grows exponentially with level. The other four strategies are total- 
order expansions of increasing order based on the following quadrature rules: linear growth 
Gauss-Legendre, exponential growth Gauss-Legendre, Clenshaw-Curtis, and Gauss-Patterson. 
All the rules were selected so the final rule would have around 10^ points. 
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We consider 20 random realizations of each Genz function in rf = 5 dimensions; random 
parameters for the Genz functions were generated as in [19]. Each estimate of LP' approximation 
error is computed by Monte Carlo sampling with 10'^ points. Figure 9.1 plots LP error at each 
stage, where each point represents the mean error over the 20 random functions. 

We can draw relatively simple conclusions from this data. All the methods show fast con- 
vergence for the first four continuous Genz functions, indicating that the internal aliasing issues 
have indeed been resolved. In contrast, one would expect direct quadrature to suffer from large 
aliasing errors for the three super-linear growth rules. Otherwise, judging the efficiency of the 
different rules is not prudent, because differences in truncation and the functions themselves 
obscure differences in efficiency. In deference to our adaptive strategy, we ultimately do not 
recommend this style of isotropic and function-independent truncation anyway. We do note, 
however, that full tensor rules never perform better than the others; this is to be expected, as 
the Genz functions are not completely coupled. The last two functions are not continuous, so all 
the methods converge much more slowly, again as expected. 

To test our adaptive approach. Figure 9.2 shows results from a similar experiment, now 
comparing the convergence of an adaptive Smolyak pseudospectral algorithm with that of a 
non-adaptive algorithm. For consistency, both are based on Gauss-Patterson quadrature. As 
we cannot synchronize the number of evaluations used by the adaptive algorithm for different 
functions, we plot individual errors for the 20 random functions instead of the mean error. This 
reveals the variability in difficulty of the functions, which was hidden in the previous plot. We 
conclude that the adaptive algorithm also converges as expected, with performance comparable 
to the non-adaptive algorithm. We do not expect much improvement from adaptivity here, 
because the Genz functions are essentially isotropic with smoothly decaying coefficients, so there 
is little structure for the adaptive algorithm to exploit. Although omitted here for brevity, other 
quadrature rules produce similar results. 

9.2. Adaptivity: chemical kinetics. To illustrate the benefits of an adaptive Smolyak 
approach, we build a surrogate for a realistic simulation of a combustion kinetics problem. Specif- 
ically, we consider the auto-ignition of a methane-air mixture given 14 uncertain rate parameters. 
Governing equations for this process are a set of stiff nonlinear ODEs expressing conservation 
of energy and of chemical species [18]. The uncertain rate parameters represent activation ener- 
gies of reactions governing the conversion of methane to methyl, each endowed with a uniform 
distribution varying over [0.8, 1.25] of the nominal value. These parameters appear in Arrhenius 
expressions for the species production rates, with the reaction pathways and their nominal rate 
parameters given by the GRIMech 3.0 mechanism [15]. The output of interest is the logarithm 
of the ignition time, which is a functional of the trajectory of the ODE system. Simulations 
were performed with the help of the TChcm software library [27], which provides convenient 
evaluations of thermodynamic properties and species production rates, along with Jacobians for 
implicit time integration. 

Chemical kinetics are an excellent testbed for adaptive approximation because, by the nature 
of detailed kinetic systems, we expect strong coupling between some inputs and weak coupling 
between others, but we cannot predict these couplings a priori. We test the effectiveness of adap- 
tive Smolyak pseudospectral methods based on the four quadrature rules discussed earlier. As 
our earlier analysis suggested that Gauss-Patterson quadrature should be most efficient, our basis 
of comparison is a non-adaptive Gauss-Patterson total-order Smolyak pseudospectral expansion. 
We ran the non-adaptive algorithm with a total order index set truncated at ti = 5 (which in- 
cludes monomial basis terms up through ^'23)' using around 40000 point evaluations and taking 
over an hour to run. We tuned the four adaptive algorithms to terminate with approximately 
the same number of evaluations. 

Figure 9.3 compares convergence of the five algorithms. The errors reported on the 
vertical axis are Monte Carlo estimates using 1000 points. Except for a small deviation at fewer 
than 20 model evaluations, all of the adaptive methods significantly outperform the non-adaptive 
method. The performance of the different quadrature rules is essentially as predicted in Section 
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7: Gauss-Patterson is the most efficient, exponential growth Gauss-Legendre and Clenshaw- 
Curtis are nearly equivalent, and linear growth Gauss-Legendre performs worse as the order of 
the polynomial approximation increases. Compared to the non-adaptive algorithm, adaptive 
Gauss-Patterson yields more than two orders of magnitude reduction in the error at the same 
number of model evaluations. Linear growth Gaussian quadrature is initially comparable to 
exponential growth Gaussian quadrature, because the asymptotic benefits of exponential growth 
do not appear while the algorithm is principally using very small one-dimensional quadrature 
rules. At the end of these experiments, a reasonable number of higher order quadrature rules are 
used and the difference becomes visible. 

We conclude by illustrating that the adaptive algorithm is effective because it successfully 
focuses its efforts on high-magnitude coefficients — that is, coefficients that make the most sig- 
nificant contributions to the function. Even though the non-adaptive expansion has around 
37,000 terms and the final adaptive Gauss-Patterson expansion only has about 28,000 terms, 
the adaptive expansion exhibits much lower error because most of the additional terms in the 
non-adaptive expansion are nearly zero. By skipping many near-zero coefficients, the adaptive 
approach is able to locate and estimate a number of higher-order terms with large magnitudes. 
Figure 9.4 depicts this pattern by plotting the difference between the numbers of included terms 
in the final adaptive Gauss-Patterson and non-adaptive expansions. The adaptive algorithm does 
not actually add any higher order monomials; neither uses one-dimensional basis terms of order 
higher than ^23- Instead, the adaptive algorithm adds mixed terms of higher total order, thus 
capturing the coupling of certain variables in more detail than the non-adaptive algorithm. The 
figure shows that terms through 30*^ order are included in the adaptive expansion, all of which 
are products of non-constant polynomials in more than one dimension. 

10. Conclusions. This paper gives a rigorous development of Smolyak pscudospcctral algo- 
rithms, a practical approach for constructing polynomial chaos expansions from point evaluations 
of a function. A common alternative approach, direct quadrature, has previously been shown 
to suffer from large errors. We explain these errors as a consequence of internal aliasing and 
delineate the exact circumstances, derived from properties of the chosen polynomial basis and 
quadrature rules, under which internal aliasing will occur. Internal aliasing is a problem inherent 
to direct quadrature approaches, which use a single (sparse) quadrature rule to compute a set of 
spectral coefficients. These approaches fail because they substitute a numerical approximation 
for only a portion of the algorithm, i.e., the evaluation of integrals, without considering the im- 
pact of this approximation on the entire construction. For almost all sparse quadrature rules, 
internal aliasing errors may be overcome only through an inefficient use of function evaluations. 
In contrast, the Smolyak pseudospectral algorithm computes spectral coefficients by assembling 
tensor-product pseudospectral approximations in a coherent fashion that avoids internal aliasing 
by construction; moreover, it has smaller external aliasing errors. To establish these properties, 
we prove that the accurate set of a Smolyak pseudospectral approximation contains a union of 
the accurate sets of all its constituent tensor-product approximation operators. These results are 
applicable to any choice of quadrature rule and generalized sparse grid, and are verified through 
numerical demonstrations. 

A key strength of Smolyak algorithms is that they are highly customizable through the choice 
of admissible multi-index sets. To this end, we describe a simple alteration to the adaptive sparse 
quadrature approach of Gerstner & Griebel [13], creating a corresponding method for adaptive 
pseudospectral approximation. Numerical experiments then evaluate the performance of different 
quadrature rules and of adaptive versus non-adaptive pseudospectral approximation. Tests of 
the adaptive method on a realistic chemical kinetics problem show multiple order-of-magnitude 
gains in accuracy over a non-adaptive approach. 

While the adaptive approach illustrated here is deliberately simple, many extensions are 
possible. For instance, as described in Section 6.3, measures of computational cost may be added 
to the local refinement criterion. One could also use the gradient of the error indicator to 
identify optimal directions in the space of multi-indices along which to continue refinement, or 
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to avoid adding all the forward neighbors of the multi-index selected for refinement. These and 
other developments will be pursued in future work. 

A flexible open-source C-h- code implementing the adaptive approximation method discussed 
in this paper is available at https://bitbucket.org/mituq/muq/. 
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Figure 9.1. Mean convergence of the non-adaptive isotropic total-order Smolyak pseudospectral algorithm 
with various quadrature rules, compared to the full tensor pseudospectral algorithm, on the Genz test functions. 
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Figure 9.2. convergence of the adaptive and non-adaptive Gauss-Patterson Smolyak pseudospectral 

algorithm. Individual results for 20 random instances of the Genz functions are shown. 



ADAPTIVE SMOLYAK PSEUDOSPECTRAL APPROXIMATIONS 



27 




Number of Evaluations 



Figure 9.3. 

convergence of ignition delay in a 14- dimensional chemical kinetic system; comparing 
a {non-adaptive isotropic total-order Gauss-Patterson-based Smolyak pseudospectral algorithm} to the adaptive 
algorithm with various quadrature rules. 




Total Order of Basis 

Figure 9.4. The plot depicts the difference betvieen the number of coefficients of a particular magnitude 
and order in the final adaptive and non-adaptive Gauss-Patterson based expansions. The horizontal axis is the 
order of the term and the vertical axis specifies logj^g "/ ^■^^ coefficient value. The color represents logj^Q of the 
difference between the two methods, where positive values indicate more terms in the non-adaptive expansion. 
Hence, the dark blue at (6,-10) indicates that the non-adaptive expansion includes around 3,000 extra terms of 
magnitude 10"^" and the dark red at (10,-8) indicates that the adaptive expansion includes about 1,000 extra 
terms of magnitude 10~*. Grey squares are the same for both expansions and white squares are not present in 
either. 



